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Abstract 

We carry out state of the art simulations of properties of confined liquid helium near the 
superfluid transition to a degree of accuracy which allows to make predictions for the outcome of 
fundamental physics experiments in microgravity. First we report our results for the finite-size 
scaling behavior of heat capacity of superfluids for cubic and parallel-plate geometry. This allows 
us to study the crossover from zero and two dimensions to three dimensions. Our calculated 
scaling functions are in good agreement with recently measured specific heat scaling functions 
for the above mentioned geometries. We also present our results of a quantum simulation of sub- 
monolayer of molecular hydrogen deposited on an ideal graphite substrate using path-integral 
quantum Monte Carlo simulation. We find that the monolayer phase diagram is rich and very 
similar to that of helium monolayer. We are able to uncover the main features of the complex 
monolayer phase diagram, such as the commensurate solid phases and the commensurate to 
incommensurate transition, in agreement with the experiments and to find some features which 
are missing from the experimental analysis. 


1 Specific Heat Scaling Function 

First, we present the results of our simulation of the specific heat scaling function for superfluids 
confined in cubic and film geometries. 

Even though earlier experiments on superfluid helium films of finite thickness [2] seemed to 
confirm the validity of the finite-size-scaling(FSS), there were later experiments[3, 4] where it was 
shown that the superfluid density of thick helium films does not satisfy FSS when the expected 
values of critical exponents were used. Similarly, in measurements of the specific heat of helium in 
finite geometries, other than the expected values for the critical exponents were found [5]. 

More recent experiments in nricrogravity environment [6] as well as earth bound experiments [7, 8] 
are consistent with scaling and they have determined the specific heat scaling function for the paral- 
lel plate (film) geometry (case (a)) and they are in reasonable agreement with the scaling function as 
was predicted by Monte Carlo simulations [9, 10] and renormalization group techniques [11]. While 
the specific heat scaling function for case (b) confinement has been theoretically determined[12] 
and it was found to be significantly suppressed relative to the case (a) there are so-far no experi- 
mental data to compare. More recently, the specific heat scaling function for the case (c) has been 
experimentally determined[13, 14]. 

The main goal of this part of this report is to present the results of our Monte Carlo simulations 
to determine the specific heat scaling function for cubes with open boundary conditions (BC) in 
all three directions (confining case (c)). In this case the scaling function characterizes the zero- 
dimensional to three-dimensional transition. Our results for the scaling function are compared to 
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the very recently obtained experimental results for specific heat scaling function in the case of 
cubic confinement [13, 14]. We find satisfactory agreement with no free parameters. In addition, we 
present results for the specific heat scaling function for the parallel plate geometry on lattices of 
size L\ x L 2 x L with L\ = L^» L where we have applied periodic BC along the £ 1 ,2-directions 
and open BC along the film-thickness direction of size L. The latter case was carried out in order to 
compare the results for Dirichlet BC (on the top and bottom of the film) obtained earlier[9, 10]. In 
Refs [9, 10] it was found that while the results with periodic BC along the film-thickness direction 
were very different from those obtained with Dirichlet BC, the results obtained with Dirichlet BC 
fit the experimental results with no free parameter. In this paper we find that the scaling function 
obtained with open BC along the finite dimension is close to that obtained with Dirichlet and 
also fits reasonably well the experimental results obtained by the Confined Helium Experiment [6] 
(CHEX). 

We have performed a numerical study of the scaling behavior of the specific heat of 4 He in a 
cubic and in a film geometry at temperatures close to the critical temperature T\. The superfluid 
transition of liquid l He belongs to the universality class of the three-dimensional x — y model, thus, 
we are going to use this model to compute the specific heat at temperatures near T\ using the cluster 
Monte-Carlo method [15]. Using Wolff’s cluster Monte Carlo updating algorithm [15], we computed 



Figure 1: The computed universal function f 2 (x) for open and periodic BC and for cubes is com- 
pared to the experimental results[13, 14]. We have used the experimental values for the critical 
exponents [17] to avoid the use of fitting parameters. 


the specific heat c(T, L ) of the x — y model as a function of temperature on several cubic lattices L 3 
(with L = 20,30,40,50 ). Open (free) boundary conditions were applied in all directions, namely 
the spins on the surface of the cube are free to take any value. These surface spins interact only with 
the 5 nearest neighbors, one in the interior and 4 on the surface of the cube and there is one missing 
neighbor. The results for the specific heat scaling function / '2 for this geometry[19] and boundary 
conditions are compared with the experimental results in Fig. 1. Taking into consideration the 
fact that we have used no fitting parameter the agreement is excellent. To further investigate the 
role of open BC, we have also calculated the specific heat scaling function f\ (x) for the case of the 
parallel plate geometry L\ x L 2 x L (L 1,2 >> L) using periodic boundary conditions along the long 
directions of the film and open BC along the thickness direction L. For this case we need to take 
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Figure 2: Film geometry: The computed universal function fi(x) with open BC (solid circles) is 
compared to the previously calculated scaling function using Dirichlet BC[9](data shown as stars) 
and periodic BC (shown as plus signs) and the experimental results of Lipa et al.[ 6] (open circles) 
and those of Mehta et al. [8] (open triangles). 


the limit L\ 2 — ► 00 first; in Ref. [10] it was found that using L\ = L -2 = 5 L was large enough, in the 
sense that systematic errors due to the finite-size effects from the fact that L\ 2 are not infinite are 
smaller than the statistical errors for realistic computational time scales. The present simulations 
for films were done on lattices 60 x 60 x 12, 70 x 70 x 14, and 80 x 80 x 16. In Fig. 2 our results for films 
with various boundary conditions are compared with the experimental results. Notice that while 
periodic boundary conditions disagree with the experimental scaling function when more realistic 
boundary conditions, such as open or Dirichlet, are used the agreement becomes satisfactory. 

2 Sub-monolayer H 2 on graphite 

In this part of our report we present our results of our quantum simulation studies [20] of sub- 
monolayer molecular hydrogen on graphite in our effort to understand quantum films (helium and 
hydrogen) . 

The commensurate-incommensurate transition has been extensively studied in classical mono- 
layer films both theoretically and experimentally. Monolayers of hydrogen or helium are quantum 
mechanical systems and, in principle, one might suspect a different behavior due to the effect of 
strong quantum fluctuations. 

The monolayer phase diagram of molecular para- hydrogen adsorbed on graphite as inferred from 
experimental studies[21] is shown in Fig. 3. This phase diagram was originally drawn based on the 
anomalies found in the specific heat [22, 23, 24]; for the characterization of the various phases, low 
energy electron diffraction (LEED)[26, 27] as well as neutron diffraction[28, 21, 29] studies have 
been carried out. At 1/3 coverage the molecules condense on the surface of graphite in a \/3 x \/3 
commensurate solid. In the low coverage region (p < 0.6) (in units of the \/3 x y/3 commensurate 
density), it forms a commensurate solid-gas coexistence phase at low temperature and a 2D gas 
phase at higher temperature (T > 10/i). For coverages 0.6^p^0.9, as a function of temperature, 
there is a transition from a commensurate solid cluster phase to a commensurate solid phase with 
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vacancies at higher temperature and to a 2D gas at even higher temperature. At density somewhat 
higher than unity the so-called a-phase is formed which is believed to be a striped domain-wall 
solid phase and at higher temperature it transforms to the so-called /3-phase and to a fluid phase at 
even higher temperature. At high densities a transition take place to a triangular incommensurate 
solid phase. We have recently studied [32] the phase diagram of molecular hydrogen on graphite at 



Figure 3: Phase diagram of molecular hydrogen adsorbed on graphite from Ref. [29] reproduced 
here for easy reference. Density 1.0 corresponds to the complete commensurate \/3 x y/S phase. 


and below 1/3 coverage using Path Integral Monte Carlo (PIMC) simulation. Our computed phase 
diagram was in general agreement with that inferred from the experimental studies for this coverage 
range. Here we report our results obtained by PIMC calculation to study the phase diagram above 
the commensurate density which is less well understood both theoretically and experimentally. 
At coverages above p ~ 1.05 the monolayer undergoes a commensurate-incommensurate(C-IC) 
transition. While the existence of the so-called a and /? phases were known from the specific 
heat anomalies, their characterization came from LEED[26, 27] and neutron scattering[21] studies. 
Using LEED experiments Cui and Fain[26, 27] observed a uniaxial IC solid with striped superheavy 
domain walls and a rotated triangular IC solid at higher coverages p > 1.23. Freimuth, Wiechert, 
and Lauter[21] (FWL) presented a neutron diffraction study of the C-IC transition and their results 
are in agreement with the LEED results. They examined the commensurate-incommensurate solid 
transition, especially the striped domain-wall solid phase. In the a-phase the diffraction intensity 
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has a main peak characteristic of a compressed lattice and a satellite peak on the lower side of the 
commensurate peak wave- vector position k c = 1.703 A -1 that arise from the spatial modulation due 
to ordered striped-domain walls. As coverage increases, the separation of the two peaks increases 
and the satellite intensity decreases. At higher temperatures the peak height drops and the satellite 
peak vanishes which implies that the commensurate domains vanish and the system becomes an 
isotropic fluid phase. As coverage increases, the first molecular hydrogen layer forms a rotated 
incommensurate solid (RIC) phase. Neutron diffraction studies show a sharp and intense peak 
at wave-vector k = l^A^ 1 at the density p = 1.34 which reveals an IC equilateral triangular 
structure. This phase is continuously compressed as the coverage increases up to the highest 
density allowed before layer promotion. 

In parallel to these experimental investigations several important theoretical studies were under- 
taken to understand the commensurate-incommensurate transition [33, 34, 35, 36, 37] in physiorbed 
systems. Molecular hydrogen and atomic helium physisorbed on graphite are quantum films char- 
acterized by strong zero point motion. Therefore, one could question the degree to which a classical 
picture might be valid and might expect new phenomena and phases to occur. Motivated by these 
thoughts we extended our earlier investigation [32] to study this system above 1/3 coverage up 
to layer completion starting from the known hydrogen molecule-molecule and hydrogen-graphite 
interactions [40, 41, 42, 43] and using the PIMC[39] which is a Quantum Monte Carlo technique 
adopted for the study of strongly interacting quantum films by Pierce and Manousakis[38]. 



(A) (B) 

Figure 4: A(B) Contour plot of the probability distribution at the density 0.0716A -2 (0.0694 A -2 ) 
and T= 1.33 K. The simulation cell contains 60 (80 molecules). The dots indicate graphite 
adsorption sites. The solid structure is uniaxially compressed along our y direction. In addition, 
in both case (A) and (B), two commensurate solid domains separated by the incommensurate solid 
domains can be seen which implies a mass density wave along the same direction. 


We used simulation cells that can accommodate the structures that have been proposed by 
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FWL. a) First, we consider p = 0.0694 A -2 . We have chosen a simulation cell with dimensions 
x= 5\/3 a gr and y= 22 a gr , where a gr = 2.459A is the carbon-carbon distance on the graphite 
surface and 80 hydrogen molecules to produce the above density. In Fig. 4(B) we give the contour 
plot of the calculated probability distribution where we see that two commensurate-solid domains 
are separated by the incommensurate solid domains. The solid structure is uniaxially compressed 
along our y direction such that a new row of molecules is introduced for every 8 rows. Notice that 
the period in the x-direction is \/3 a gr and there is modulation along the y direction with period 
llcig r , so that the wavelength A s of this striped domain- wall modulation is A s = 27.049. b) Second, 
we have performed the simulation at the density p = 0.0716A -2 . We have chosen a simulation cell 
with dimensions x = 5\Z3a gr and y = 16 a gr and 60 hydrogen molecules in order to achieve the 
above mentioned density. The calculated contour plot of the probability distribution is shown in 
Fig. 4(A). Notice that in this case also there is an ordered striped domain- wall solid phase along 
our y direction. Namely, the amplitude of the molecular density wave is modulated along the y axis 
with a period of about half our cell size along the y axis. Notice that there are two commensurate 
solid domains separated by denser regions. 



Figure 5: Static structure factor S(k) at the density 0.0716JW 2 (0.0694A -2 ) for T= 1.33 K as a 
function of the magnitude of k. The main Bragg peaks are at k\= 1.759 A -1 and k 2 = 1.916 A ^ 1 
( ki= 1.743 A ' 1 and k 2 = 1.858 A -1 ). There is a satellite peak located at k c sat = 1.597 A ” 1 ( k c sat = 
1.626 A -1 .) Notice that the satellite peak at in Fig. (A) is significantly weaker than the satellite 
peak at the lower density shown in Fig. (B) 

Our computed static structure factor S(k) for p = 0.0694A~ 2 is shown in Fig. 4(A). The main 
Bragg peaks of this structure occur at ki = (1.475 A” 1 , 0.929 A -1 ) and k 2 = (0,1.858 A -1 ). The 
satellite peak occurs at k^ at = (0,1.626 A” 1 ). The experimental values of the magnitude of the 
wave-vectors at the peaks are k^ in = 1.743 A -1 and k e s ^= 1.632 A” 1 , which compare well with 
our computed values of k\= 1.743 A _1 and k c sat = 1.626 A -1 . We believe that FWL could not 
observe the peak at k 2 because of the strong interference with the (002) graphite peak. 

The interpretation of these results is as follows: Analyzing the contour plot of Fig. 4(B) we 
find that a) the solid appears uniaxially compressed along the y direction, by adding another row 
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for every 10 rows of molecules and spreading them evenly, while along the x direction the average 
spacing between the molecules remains that of the commensurate solid, b) superimposed there is 
a density modulation along the y direction which has wavelength several times greater than the 
average nearest neighbor distance but of small amplitude. The actual size of the wavelength of 
this striped domain- wall modulation is half of the length of our simulation cell in the y-direction, i.e., 
A s = 27.05^, as discussed in the previous paragraph. The two main Bragg peaks correspond to the 
unit vectors which span the reciprocal space of this uniaxially compressed triangular solid structure. 
The satellite peaks are due to lateral modulation in our y direction and are located at k sat = ki^— k s . 
k ] .2 correspond to the wave vector of each of the two main Bragg peaks and k s = (0, 27t/A s ), where 
A s is the wave length of a unit cell for the striped domain-wall solid structure in the modulated 
direction. Using this A s value, k 2 = (0, 1.858A -1 ) and k s = (0,0.232A -1 ) (obtained by using the 
values of A s = 27.049^4 found by inspection of Fig. 4(B)) we can find that the satellite peak position 
is at k sa t = (0, 1.626A -1 ), which agrees well with our peak value. 

2.1 Domain Wall Melting 

In Fig. 6 we present the contour plot of the the probability distribution where a striped domain- 
wall fluid phase (at T = 11.11 K) and a fluid phase (at T = 16.67) are evident. Notice that the 
stripe domain-wall fluid phase (left part of Fig. 6) is characterized by mobile commensurate and 
incommensurate domains. To understand this contour plot we need to be reminded of the periodic 
boundary conditions used in our simulation and the fact that these domains cannot intersect each 
other. This implies that one domain will oscillate back and forth between the neighboring domains. 

We would like to call the reader’s attention on the presence of strong finite-size effects in this 
study of the domain wall melting. In addition, the domain boundary melting is analogous to the 
roughening of vicinal surfaces and meandering of steps, and our simulation cell is not wide enough 
to encompass formation of kink-pairs along the domain boundaries. We can very approximately 
determine the melting temperature of the domain- wall solid phase from the temperature dependence 
of the specific heat and from the temperature dependence of the peak height of the static structure 
factor S(k). 



Figure 6: Contour plot of the probability distribution at the density 0.0716A 2 at two temperatures 
11.11 K and 16.67 K. 


The calculated specific heat as a function of temperature is shown in Fig. 7 and is characterized 
by two peaks, the first corresponds to the melting of the domain- wall solid while the second peak 
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0 10 20 30 

Temperature (K) 

Figure 7: Specific heat versus temperature at the density O.OTIG^ -2 . The long-dashed line is a spline 
fit to the specific heat values and is a guide to the eye. The filled circles are the experimental specific 
heat at the density 0.0716A -2 . The simulation cell is 29.813 A x 39.344 A (84 H2 molecules). Our 
computed values for the melting and evaporation temperature are the peak positions, 9.09 K and 
12.5 K , respectively. 
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indicates the melting of the stripes and the solid into a fluid. Our computed values for the melting 
temperature are generally lower that the experimental values. For example, for p = 0.0716A~ 2 
we find that at T = 9.09/i the striped domain-wall solid phase undergoes a transition to the 
domain-wall fluid phase and at T = 12.5 K the /3-phase becomes an isotropic fluid phase. Factors 
for obtaining lower values for the critical temperature than the experimental values could be the 
finite-size effects and the interaction strength used to describe the hydrogen-graphite interaction. 
We also computed the static structure factor shown in Fig. 8 and we studied the temperature 
dependence of the height of its first main peak. The peak height significantly decreases near the 
melting temperature of the domain-wall solid and near the domain-wall evaporation temperature. 

Clearly much larger size lattices are needed to study the nature of this phase transition. Such 
a study is beyond our current computational resources and we restrict ourselves to our qualitative 
results obtained here with full quantum mechanical treatment of the problem. In a different line 
of work one can use classical simulations of models of such stripe melting to carefully examine 
finite-size effects. 

2.2 Rotated Incommensurate Solid 

The first layer of molecular hydrogen adsorbed on graphite forms an incommensurate solid phase 
before the first layer is complete. We have carried out a simulation at the density p = 0.0849Ji 2 
using a simulation cell with dimensions x= 9 \/3 a gr and y= 9 a gr . The calculated probability 
distribution(Fig. 9) also clearly shows the presence of an equilateral triangular solid structure 
which is not registered with the underlying graphite lattice and it is rotated relative to the graphite 
lattice. This rotation was first predicted and discussed by Novaco and McTague[33]. The angle of 
the incommensurate lattice relative to the graphite lattice is approximately 5° in good agreement 
with the experimental value[27] for the above density. The calculated static structure factor has 
a single sharp peak at k= 1.99 Ar 1 corresponding to an unregistered equilateral triangular lattice 
in agreement with the value of k= 1.97 A reported from the neutron diffraction experimental 
study [21], 



Figure 9: Contour plot of the probability distribution at the density 0.0849A 2 . The simulation 
cell is 38.331 A x 22.131 A (72 H2 molecules). 
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